Dynamics of stick-slip in peeling of an adhesive tape 
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We investigate the dynamics of peeling of an adhesive tape subjected to a constant pull speed. 
We derive the equations of motion for the angular speed of the roller tape, the peel angle and the 
pull force used in earlier investigations using a Lagrangian. Due to the constraint between the 
pull force, peel angle and the peel force, it falls into the category of differential-algebraic equations 
requiring an appropriate algorithm for its numerical solution. Using such a scheme, we show that 
stick-slip jumps emerge in a purely dynamical manner. Our detailed numerical study shows that 
these set of equations exhibit rich dynamics hitherto not reported. In particular, our analysis 
shows that inertia has considerable influence on the nature of the dynamics. Following studies in 
the Portevin-Le Chatelier effect, we suggest a phenomenological peel force function which includes 
the influence of the pull speed. This reproduces the decreasing nature of the rupture force with 
the pull speed observed in experiments. This rich dynamics is made transparent by using a set 
of approximations valid in different regimes of the parameter space. The approximate solutions 
capture major features of the exact numerical solutions and also produce reasonably accurate values 
for the various quantities of interest. 

PACS numbers: 05.45.Pq, 62.20.Mk, 68.35.Np 



I. INTRODUCTION 



Peeling is a kind of fracture that has been stud- 
ied experimentally in the context of adhesion and is a 
technologically important subject. Experimental stud- 
ies on peeling of an adhesive tape mounted on a cylin- 
drical roll are usually in constant pull speed condition 
IB IE [3 IE El- More recently constant load experi- 
ments have also been reported [E ■ Early studies by 
Bikermann 5] , Kaeble |6| have attempted to explain the 
results by considering the system as a fully clastic object. 
This is clearly inadequate as it ignores the viscoelastic 
nature of the glue at the contact surface and therefore 
cannot capture many important features of the dynamics. 
The first detailed experimental study of Maugis and Bar- 
quins show stick-slip oscillations within a window of 
pull velocity with decreasing amplitude of the pull force 
as a function of the pull velocity. Further, these authors 
report that the pull force shows sinusoidal, sawtooth and 
highly irregular (chaotic as these authors refer to) wave 
patterns with increasing velocities. More recently, Gan- 
dur et al. have carried out a dynamical time series anal- 
ysis of the force waveforms, as well as those of acoustic 
emission signals and report chaotic force waveforms at 
the upper end of the pull velocities Q • One characteris- 
tic feature of the peeling process is that the experimental 
strain energy release rate shows two stable branches sep- 
arated by an unstable branch. Stick-slip behavior is com- 
monly observed in a number of systems such as jerky flow 
or the Portevin-Le Chatelier (PLC) effect Q, frictional 
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sliding 0, and even earthquake dynamics is thought to 
result from stick-slip of tectonic plates ^(j- Stick-slip 
is characterized by the system spending most part of 
the time in the stuck state and a short time in the slip 
state, and is usually seen in systems subjected to a con- 
stant response where the force developed in the system is 
measured by dynamically coupling the system to a mea- 
suring device. One common feature of such systems is 
that the force exhibits "negative flow rate characteristic" 
(NFRC). Models which attempt to explain the dynamics 
of such systems use the macroscopic phenomenological 
NFRC feature as an input, although the unstable region 
is not accessible. This is true for models dealing with the 
dynamics of the adhesive tape as well. To the best of our 
knowledge, there is no microscopic theory which predicts 
the origin of the NFRC macroscopic law except in the 
case of the PLC effect [11 E2 (see below). 

As there is a considerable similarity between the peel- 
ing of an adhesive tape and the PLC effect, it is useful to 
consider the similarities in some detail. The PLC effect 
refers to a type of plastic instability observed when sam- 
ples of dilute alloys are deformed under constant cross 
head speeds [l3j . The effect manifests itself in the form 
of a series of serrations in a range of applied strain rates 
and temperatures. This feature is much like the peeling 
of an adhesive tape. Other features common to these 
two situations are: abrupt onset of the large amplitude 
oscillations at low applied velocities with a gradually de- 
creasing trend and NFRC, which in the PLC effect refers 
to the existence of negative strain rate sensitivity of the 
flow stress. In the case of the PLC effect, the physical 
origin of the negative strain rate sensitivity is attributed 
to the ageing of dislocations and their tearing away from 
the cloud of solute atoms. Recently, the origin of the 
negative SRS has been explicitly demonstrated as aris- 
ing from competing time scales of pinning and unpinning 
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in the Ananthakrishna's model [lit . In the case of 
adhesive tape, the origin of NFRC can be attributed to 
the viscoelastic behavior of the fluid. (Constant load and 
constant load rate experiments are possible in the PLC 
also.) While simple phenomenological models based on 
NFRC explain the generic features of the PLC effect > 
there appears to be some doubts if the equations of mo- 
tion conventionally used in the present case ofpeeling are 
adequate to describe the velocity jumps 0, E|- Indeed, 
these equations of motion are singular and pose problems 
in the numerical solutions. 

Apart from detailed experimental investigation of the 
peeling process, Maugis and Barquins 0, have also con- 
tributed substantially to the understanding of the dy- 
namics of the peeling process. However,the first dy- 
namical analysis is due to Hong and Yue 0] who use an 
"N" shaped function to mimic the dependence of the peel 
force on the rupture speed. They showed that the sys- 
tem of equations exhibits periodic and chaotic stick-slip 
oscillations. However, the jumps in the rupture speed 
are introduced externally once the rupture velocity ex- 
ceeds the limit of stability pi Il5|. Thus, the stick-slip 
oscillations are not obtained as a natural consequence of 
the equations of motion. Therefore, in our opinion the 
results presented in Ref. Q are the artifacts of the nu- 
merical procedure followed. Ciccotti et al. interpret 
the stick-slip jumps as catastrophes. Again, the belief 
that the jumps in the rupture velocity cannot be ob- 
tained from the equations of motion appears to be the 
motivation for introducing the action of discrete opera- 
tors on the state of the system to interpret the stick-slip 
jumps though they do not demonstrate the correct- 
ness of such a framework for the set of equations. Lastly, 
there are no reports that explain the decrease in the am- 
plitude of the peel force with increasing pull speed as 
observed in experiments. As there is a general consen- 
sus that these equations of motion correctly describe the 
experimental system, a proper resolution of this question 
( on the absence of dynamical jumps in these equations ) 
assumes importance. 

The purpose of this paper is to show that the dynam- 
ics of stick-slip during peeling can be explained using a 
differential-algebraic scheme meant for such singular sit- 
uations and demonstrate the rich dynamics inherent 
to these equations. In what follows we first derive the 
equations of motion (used earlier 0) by introducing an 
appropriate Lagrangian for the system. Then, we use 
an al gori thm meant to solve differential-algebraic equa- 
tions [l6| and present the results of our simulations for 
various parameter values. One of our major findings is 
that inertia has a strong influence on the dynamics. In 
addition, following the dynamization scheme similar to 
the one used in the context of the PLC effect 0], we 
suggest that the peel force depends on the applied veloc- 
ity. Using this form of peel force leads to the decreasing 
nature of the magnitude of the pull force as a function 
of applied velocity. For certain values of the inertia, we 
find canard type solutions. These numerical results are 




FIG. 1: Schematic plot of experimental setup 



captured to a reasonable accuracy using a set of approxi- 
mations valid in different regimes of the parameter space. 
Even though, our emphasis is on demonstrating the cor- 
rectness of these equations of motion and richness of the 
inherent dynamics that capture the qualitative features 
of the peeling process, we also attempt to make a com- 
parison of the experimental results mentioned above to 
the extent possible. 

II. EQUATIONS OF MOTION 

For the sake of completeness, we start by considering 
the geometry of the experimental setup shown schemat- 
ically in Fig. ^ An adhesive roll of radius R is mounted 
on an axis passing through O normal to the paper and 
is pulled at a constant velocity V by a motor positioned 
at O' with a force F acting along PO'. Let the distance 
between O and O' be /, and that between the contact 
point P to O' be L. The point P moves with a local 
velocity v which can undergo rapid bursts in the velocity 
during rupture. The force required to peel the tape is 
usually called the force of adhesion denoted by /. The 
two measured branches referred to earlier, are those of 
the function / in a steady state situation of constant 
pulling velocity (i.e., there are no accelerations). The 
line L makes an angle 9 with the tangent at the contact 
point P. The point P subtends an angle a at O, with 
the horizontal line OO' . We denote the elastic constant 
of the adhesive tape by fc, the elastic displacement of the 
tape by u, the angular velocity by uj and the moment 
of inertia of the roll by /. The angular velocity itself is 
identified by lu = a + v / R. The geometry of the setup 
gives L cos 9 = —I sin a and L sin9 — I cos a — R which 
further gives, L 2 = l 2 + R 2 — 21R cos a. The total velocity 
V at O' is then made up of three contributions , given 
by V = v + u — L, which gives 

v = V + L - ii = V - R cos9 a - u. (1) 

Following standard methods in mechanics, it is straight- 
forward to derive the equations of motion for a and lo 
by considering (a, d, u, ii) as the generalized co-ordinates. 
The corresponding Lagrangian of the system can be writ- 
ten as 

I fa 
C(a,a,u,u) = ~[uj(a,d,u,u)} 2 - -u 2 . (2) 
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We write the dissipation function as 



K = *(v,V) = / f{v,V)dv, 



(3) 



where f(v, V) physically represents the peel force which 
we assume is dependent on rupture speed as well as the 
pull speed assumed to be derivable from a potential func- 
tion $(u,U). The physical origin of this is due to the 
competition between the internal relaxation time scale of 
the viscoelastic fluid and the time scale determined by 
the applied velocity 0. When the applied velocity is 
low, there is sufficient time for the viscoelastic fluid to 
relax. As we increase the applied velocity, the relaxation 
of the fluid gets increasingly difficult and thus behaves 
much like an elastic substance. The effect of compe ting 
time scales is well represented by Deborah number |l8j 
which is the ratio of time scale for structural relaxation 
to the characteristic time scale for deformation. Indeed, 
in the studies on Hele-Shaw cell with mud as the viscous 
fluid, one observes a transition from viscous fingering to 
viscoelastic fracturing [19j with increasing rate of inva- 
sion of the displacing fluid. 

As stated in the Introduction, the existing models do 
not explain the decreasing amplitude of pull force. Simi- 
lar feature observed in the PLC serrations has been mod- 
eled using a scheme referred to as dynamization of the 
negative strain rate sensitivity (SRS) of the flow stress 
f(ep) 0,113, where e p is the plastic strain rate. Based 
on arguments similar to the preceding paragraph, they 
modify this function to depend on the applied strain rate, 
e a , i.e., the negative SRS of the flow stress is taken to be 
f(ep,e a ) such that the gap between the maximum and 
the minimum of the function f(e p , e a ) decreases with in- 
creasing i a . Following this, we consider / to depend on 
V also, in a way that the gap in / decreases as a function 
of the pull speed V (Fig. 0). 

Using the Lagrange equations of motion, 
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we obtain the same set of ordinary differential equations 
as in Ref. Q given by 

a = ui — v/R, (6) 

Iuj = FR cos 9 = -FRsina ~ -FRa (7) 

F = ku = k(V-v) - k cos 6(uR-v), (8) 

~ k[V-v + Raa], (9) 

with an algebraic constraint 

F(l - cos 6) - f(v, V) ~ F(l + a) - f(v, V) = 0. 



(10) 
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FIG. 2: (a) Plots of f(v, V) as a function of v (x axis in log 
scale) for V = 1 (solid curve), V — 2 (dashed curve), V = 4 
(dashed and dotted curve), V — 6 (dotted curve); see Eq. 
1141 . (b) Experimental strain energy release rate, G(v) curve 
as in Ref. [Units of f(v, V) is in N, G(v) in J/m 2 and v, 
V are in m/s ] 



we have used cos9 ~ —sina ~ —a. While Eqs. ©- 
© are differential equations, Eq. I|1(J|) is an algebraic 
constraint necessitating the use of differential-algebraic 
scheme to obtain the numerical solution [l6| . 

The fixed point of Eqs. ©, 0, ©, and (JTUJl is given 
bya = 0,w = V/R, v = V, F = f(V, V). (For numerical 
solution, in the above equations we have actually used 
sina in place of a.) This point is stable for f'(V, V) > 
and unstable for f'(V, V) < 0. As V is varied such that 
the sign of f'(V, V) changes from negative to positive 
value, the system undergoes a Hopf bifurcation and a 
limit cycle appears. The limit cycles reflect the abrupt 
jumps between the two positive slope branches of the 
function f(v,V). 



III. ALGORITHM 

The singular nature of these equations becomes clear 
if one were to consider the differential form of Eq. I|1(J|) 
given by 



1 



f'(v,V) 
[F(l + a)+Fa}/f, 



F(l-cos0) + Fsin( 



(11) 
(12) 



(The last equation results from the elimination of two 
second order equations for a.) In Eqs. Q, ©, and i|l(J|) 



where the prime denotes the derivative with respect to 
v. Equations (TTJ witn E q- ©, 0, and © [or ©] 
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constitute the full set of evolution equations for the vector 
(a,u),F, v). However, it is clearly singular at points of 
externum of f(v, V), requiring an appropriate numerical 
algorithm. 

We note that Eqs.®, 0, ©, and (HJ can be written 

as 
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where X = (a, lu, F,v), is a vector function that governs 
the evolution of X and M is a singular "mass matrix 
given by, 
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Equation (|13fl is a differential-algebraic equation (DAE) 
and can be solved using the so called singular pertur- 
bation technique |l(lj | in which the singular matrix M is 
perturbed by adding a small constant e such that the 
singularity is removed. The resulting equations can then 
be solved numerically and the limit solution obtained as 
e — > 0. We have checked the numerical solutions for e 
values ranging from 10~ 7 to 10~ 15 in some cases and the 
results do not depend on the value of e used as long as 
it is small. The results presented below, however, are for 
e = 10~ 7 . We have solved Eq. (|T3|> using a standard 
variable-order solver, MATLAB ODE15S program. 
We have parametrized the form of f(v, V) as 

f(v,V) = 400w - 35 + 110v 015 + lSOe^/ 11 ) - 2V 1 - 5 

-(415 - 45V QA - 0.35V~ 21f > a5 , (14) 

to give values of the extremum of the peel velocity that 
mimic the general form of the experimental curves 0. 
The measured strain energy release rate G(V) from sta- 
tionary state measurements is shown in Fig. 2(b). The 
decreasing nature of the gap between the maximum and 
minimum of f(y, V) for increasing V is clear from Fig. 
2(a). [The values of f(v, V) could not be correctly deter- 
mined as G(V) is in J/m 2 requiring more details. How- 
ever, the value of F max is closer to Ref. and the jumps 
in v are similar to those in experiments.] The reason for 
using the form given by Eq. I|14|l is that the effects of dy- 
namization are easily included through its dependence on 
the pulling velocity while more complicated terms are re- 
quired to mimic completely the experimental curve (par- 
ticularly the flat portion). However, we stress that the 
trend of the results remains unaffected when the actual 
experimental curve is used except for the magnitude of 
velocity jumps and the force values. 



IV. RESULTS 

We have studied the dynamics of the system of 
equations for a wide range of values of the parameters. 




FIG. 3: (a) A typical phase space trajectory in the v — F 
plane for V = 0.4,1 = 10 -5 . The corresponding f(v,V) is 
shown by a solid curve, (b) A phase space trajectory in the 
v — F plane for V = 1.0 and I = 1(T 5 . (c) A plot of a(t) for 
V = 1 and I = 10" 5 . (d) A plot of F(t) (period 4) for V = 2 
and / = 1CP 5 . (Units of v, V are in m/s, F in N, I in kg m 2 
and t in s.) 



We have found that transients for some regions of 
parameters space take considerable time to die out. 
The results reported here are obtained after these long 
transients are omitted. These equations exhibit rich 
dynamics, some even unanticipated. Here we report 
typical results for two important parameters, namely 
the pull velocity V (m/s) and the inertia I (kg m 2 ), 
keeping the elastic constant of the tape k — 1000 N/m, 
R = 0.1 m and I = 1 m 0. The influence of k will also 
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be mentioned briefly. ( Henceforth, we drop the units for 
the sake of brevity.) We find that the observed jumps 
of the orbit in the v-F plane occur in a fully dynamical 
way. More importantly, we find all the three possibilities 
namely, the orbit can jump when it approaches the 
limit of stability, before or beyond that permitted by 
f(v,V). The dynamics can be broadly classified into 
low, intermediate and high regimes of inertia. 

(i) Low inertia. Here also, there are three regimes: 
low, intermediate, and high pull velocity. 

(a) Consider keeping inertia I at a low value (say 
I = 10 -5 ) and V also at a low value (say, near the top, 
say V = 0.4 ). Here we observe regular saw tooth form 
for the pull force F. The phase plot in the F-v plane 
is as shown in Fig. 3(a). The corresponding function 
f(v, V) is also shown by the continuous curve. We see 
that the trajectory jumps almost instantaneously from 
B to C on reaching the maximum of f(v, V) (or from D 
to A when it reaches the minimum) . The system spends 
considerably more time on AB compared to that on CD. 
However, this feature of jumping of the trajectory at the 
limit of stability is only true for small values of I and 
when V is near the limit of stability. At slightly higher 
pull velocity, say V = 1, even for small I, say I = 10~ 5 , 
the jumps occur even before reaching the top or bottom 
( the points B and D) as can be seen from Fig. Efb) for 

V = 1. The small amplitude high frequency oscillations 
seen in the phase plots [ Fig. OJa), and Efb)] on the 
branch AB are due to the inertial effect, i.e., finite value 
of /. These oscillations are better seen on the a(t) plot 
shown in Fig. 3(c). For these values of parameters, the 
system is aperiodic. 

b) As we increase V, even as the saw tooth form of 
F is retained, various types of periodic orbits [period 4 
shown in Fig. Ofd) for V = 2] as well as irregular orbits 
are seen. In both cases (periodic as well as chaotic) the 
trajectory jumps from high velocity branch (CD) to the 
low velocity branch before traversing the entire branch 
or sometimes going beyond the values permitted by /. 
The value of F at which the orbit jumps is different 
for different cycles. For I = 10~ 5 , at high velocity, say 

V = 4, the phase plot is periodic. 

(ii) Intermediate and high inertia. 

(a) As the results of small V for intermediate and high 
inertia are similar, we illustrate the results for I = 10 -2 
and V = 1. The v — F phase plot, a, F and v are 
shown in Fig. Eta)-QJd). Consider, Fig. IHa) showing 
a typical phase space trajectory for a single cycle. The 
corresponding function f(v,V) is also shown by the 
thick continuous curve. We see that the maximum (and 
minimum) value of F is larger (or smaller) than that 
allowed by f{v,V). [This feature holds when the inertia 
is in the intermediate regime also, though the values of 
maxima (minima) of F are not significantly larger (less) 
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FIG. 4: (a) Phase space trajectory in the v — F plane for 
a single cycle for I — 10~ 2 and V = 1. The corresponding 
f(v, V) is shown by a thick solid curve, (b) Corresponding 
plots of a(t), (c) the pull force F(t) (period 8) and (d) the 
peel velocity v(t). (Units of v, V are in m/s, F in N, I in kg 
m 2 and t in s.) 



than f max (fmin)-] When the trajectory jumps from 
AB to CD at the highest value of F for the cycle, the 
trajectory stays on CD for a significantly shorter time 
compared to the small inertia case (/ = 10 ) and jumps 
back to AB well before F has reached the minimum of 
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f(v, V), i.e., AF is much smaller than f max - f min . The 
pull force F cascades down through a series of back and 
forth jumps between the two branches till the lowest 
value of F for the cycle is reached. Note that F at the 
point n is less than / m j n . For the sake of clarity, two 
different portions of the trajectory are marked abcdefg 
and ijklmna corresponding to the top and bottom 
regions of the plot. The corresponding points are also 
identified on the F(t) plot. After reaching n, the orbit 
jumps to a on AB, the trajectory decides to move up 
all the way till F reaches a maximum value (larger than 
fmax, the point 6) without jumping to the CD branch. 
This part of F as a function of time, which is nearly 
linear on AB, (i.e., the segment ab) displays a noticeable 
sinusoidal modulation. The sinusoidal form is better 
seen in a [Fig. @£b)]. Note that the successive drops in 
F are of increasing magnitude. The jumps between the 
two branches in the v — F plane are seen as bursts of v 
[Fig. EJd)]. For these values of parameters, the system 
is periodic. 

(b) As we increase V, the sinusoidal nature of F and a 
becomes more clear with its range becoming larger reach- 
ing a nearly sinusoidal at V = 4 for large I. [The range 
ab in Fig. 0Jc) expands. Compare Fig. EJa).] The mag- 
nitude of AF on the CD branch for small V and mod- 
erately or large /, gradually decreases with increasing V. 
The magnitude of AF itself decreases as I is increased. 
In the limit of large V and /, the drops in F and a be- 
come quite small which are now located near the maxima 
and minima of these curves. This is shown in Figs. Eta) 
andEIb). The sinusoidal nature is now obvious even in 
F(t) unlike for smaller V and I where it is clear only in 
a(t) for the low v branch. Note that for V = 4, the na- 
ture of f(v, 4) is nearly flat. This induces certain changes 
in the v — F phase plot that are not apparent in F and 
a. The jumps between the two branches are now concen- 
trated in a dense band at low and high values of F. In 
this case, the maximum (minimum) value of F is signifi- 
cantly larger (less) than f max (/ mm ). These rapid jumps 
between the branches manifest as jitter at the top and 
bottom of F and a. 

Unlike for small V [Fig. 01 a)], the nature of the tra- 
jectory in Fig. EI C ) is different. After reaching a critical 
value of F near the maximum value of F (the point b) , 
the orbit spirals upwards and then descends down till 
another critical value of F (the point c) is reached. Hav- 
ing reached c, the orbit monotonically comes down till d 
where it jumps to the AB branch. Beyond this point, it 
again spirals upwards till the point a is reached. There- 
after, F monotonically increases till b is reached. The re- 
gions ab and cd are the regions where F shows a near sinu- 
soidal form. The regions be and da are the regions where 
the orbit jumps between the branches rapidly. These 
manifest themselves as bursts of v which tend to bunch 
together almost into a band. [Compare Fig. E{d) with 
Fig. Efd).] It is interesting to note that the jumps be- 
tween the two branches occur exactly at points where 
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FIG. 5: (a) A plot of F{t) for V = 4 and I = 10~ 2 . (b) 
Corresponding plot of a(t). The inset shows an expanded 
plot of decreasing trend of a(t). (c) Corresponding plots of 
phase space trajectory that reflects the chaotic nature and (d) 
the peel velocity v(t). (v, V are in m/s, F in N, / in kg m 2 
and t in s.) 
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df/dv = 0, even when the maximum (minimum) of F 
are higher (lower) than that allowed by the stationary 
curve f(v,V). The variables are aperiodic for the set 
of parameters. The phase plots appear to be generated 
by an effective f(v, V) that is being cycled. [This visual 
feeling is mainly due to the fact that jumps between the 
branches still occur at the maximum and minimum of 
f(v, V).] 

The influence of k is generally to increase the range 
of the pull force F as can be easily anticipated and to 
decrease the associated time scale. 

It may be desirable to comment on the similarity of 
the nature of the force waveforms displayed by the model 
equations with those seen in experiments. As mentioned 
in the introduction, apart from qualitative statements on 
the waveforms in Ref. pj (such as periodic, sawtooth etc., 
which are seen in the model as well) , it should be stressed 
that there is a paucity of quantitative characterization of 
the waveforms. In this respect, the study by Gandur 
et al. fills the gap to some extent. These authors 
have carried out a dynamical analysis of the time series 
for various values of the pull velocities (for a fixed value 
of the inertia corresponding to their experimental roller 
tape geometry). In order to compare this result, we have 
calculated the largest Lyapunov exponent for a range of 
values of / and V. The region of chaos is in the domain 
of small pull velocities V when / is small. The maximum 
Lyapunov exponent turns out to be rather high, typically 
around 7.5 bits/s in contrast to the small values reported 
in Ref. @|. The large magnitude of the positive expo- 
nent in our case can be traced to the large changes in 
the Jacobian, as df(v, V)/dv varies over several order of 
magnitude^ 10 6 ) as a function of the peeling velocity 
and hence as a function of time. In contrast, Hong et 
al. use an N shaped curve where df(v, V)/dv is constant 
(and small) on both low and high V branches. However, 
these large values of Lyapunov exponents are consistent 
with rather high values obtained by Gandur et al. Q 
from time series analysis of the pull force. We also find 
chaos for intermediate and high inertia in the region of 
high velocities where the value of the Lyapunov expo- 
nent is small, typically 0.5. The small value here again 
can be traced to the small changes in df(v, V)/dv at high 
velocities. 

It must be mentioned that comparison with experi- 
ments is further complicated due to the presence of a two 
parameter family of solutions strongly dependent on both 
I and V. Thus, the phase diagram is complicated, i.e., 
the sequence of solutions encountered in the I-V plane 
as we change V or I or both does not in general display 
any specific ordering of periodic and chaotic trajectories 
(see Fig. 1 of Ref. [21| ) usually found in the well known 
routes to chaos. (For instance 2 n periods should be ob- 
served before the odd periods H^.) Indeed, in our model, 
we find the odd periods 3,5,7 etc, on increasing V (or J), 
without seeing all the 2™ periods. (These odd periods 
also imply chaos at parameter values prior to that cor- 
responding to these periods.) In view of this, a correct 



120 



80 



40 



O 



o 



o. 



1 



2 3 
pull velocity V 



FIG. 6: The plot shows the mean force drop AF as a function 
of the pull speed V, for two distinct values of /. The dashed 
line corresponds to / = 10~ 2 while the dotted line corresponds 
to I — 10 -5 . (v, V axe in m/s, F in N, I in kg m 2 and t in s.) 
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FIG. 7: A phase plot of canard type of solution in v — F plane 
for V = 0.4 and I = 1(T 3 . (v, V are in m/s, F in N, I in kg 
m 2 and t in s.) 



comparison with experiments requires an appropriate cut 
in the I — V plane consistent with the experimental val- 
ues of I and V even where they are given. However, as 
the values of I are not provided, full mapping of chaotic 
solutions is not possible. (We also note that Gandur et 
al. use a different tape from that used in Ref. pj, 
as is clear from the instability range, leading additional 
difficulties in comparison.) 

One quantitative result that can be compared with ex- 
periment is the decreasing trend of the force drop mag- 
nitude. We have calculated the magnitude of the force 
drops during stick-slip phase as a function the pull veloc- 
ity V for both low ( I = 1CT 5 ) and high (I = 1(T 2 ) in- 
ertia cases. Figu re |S| sho ws the monotonically decreasing 
trend of average AF(t) as V is increased, for both small 
and large /, a feature observed in experiments Q . These 
two distinct behaviors are a result of the dynamization 
of f(v,V) as in Eq. 

Finally, as another illustration of the richness of the 
dynamics seen in our numerical simulations, we show in 
Fig. [7| a plot of an orbit that sticks to unstable part 
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of the manifold before jumping back to the AB branch. 
Such solutions are known as canards [23| . Though canard 
type of solutions are rare, we have observed them for 
high values of I and low values of V. In our case, such 
solutions are due to the competition of time scale due 
to inertia and that due to v. This again illustrates the 
influence of inertia of the roll on the dynamics of peeling. 

It is clear that these equations exhibit rich and com- 
plex dynamics. A few of these features are easily un- 
derstandable, but others are not. For instance, the saw- 
tooth form of F for low inertia and low pull velocity can 
be explained as resulting from the trajectory sticking to 
stable part of f(v, V) and jumping only when it reaches 
the limit of stability. For these parameter values, as the 
time spent by the system is negligible during the jumps 
between the branches AB and CD (and vice versa), the 
system spends most of the time on the branch AB and 
much less on CD due to its steep nature. Then, from 
Eq. 0, it is clear that we should find a sawtooth form 
whenever the peel velocity v jumps across the branch to 
a value of v larger than the pull velocity V . 

However, several features exhibited by these system of 
equations are much too complicated to understand. We 
first list the issues that need to be explained. 

(I) Small I. 

(a) We find high frequency tiny oscillations superposed 
on the linearly increasing F [on the AB branch or 
better seen in the a plot Fig. E{c)]. This needs to be 
understood. 

(b) The numerical solutions show that the influence of 
inertia can be important even for small I and small V . 
For instance, the jumps between AB and CD branches 
occur even before F reaches the extremum values of /. 

(II) For intermediate and high values of inertia, for low 
V case. 

(a) We observe several relatively small amplitude saw 
tooth form of F on the descending part of the pull force 
F. These appear as a sequence of jumps between the two 
branches in the v — F plane which we shall refer to as the 
"jumping mode" . A proper estimate of the magnitude of 
AF is desirable. 

(b) In addition, there appears to be a critical value 
of F for a given cycle below which the return jumps 
from AB to CD stop and one observes a monotonically 
increasing trend in F [ab in Fig. Etc)]. 

(III) High I and high V. 

(a) The jumps between the branches occur at a very 
high frequency [Fig. EJc)] and now are located near the 
extremum values of F and a. But these regions are sep- 
arated by a stretch where the orbit monotonically in- 
creases on the AB branch and monotonically decreases 
on the CD branch. We need to elucidate the underly- 
ing causes leading to the switching between the jumping 
mode and monotonically increasing or decreasing mode. 

(b) For large V, say V = 4 and large I (Fig. [5J, the 



extent of values of F(t) range between 185 and 450 much 
beyond f(v, 4) whose range is around 300. This feature 
is less dominant for small I and small V case. 



V. APPROXIMATE ANALYSIS OF THE 
DYNAMICS 

As the dynamics is described by a coupled set of dif- 
ferential equations with an algebraic constraint, the re- 
sults are not transparent. We first attempt to get insight 
into the complex dynamics through some simple approx- 
imations valid in each of the regimes of the parameters. 
Solution of these approximate equations will require ap- 
propriate initial values for the relevant variables which 
will be provided from the exact numerical solutions. Due 
to the nature of approximations, the results are expected 
to capture only the trend and order of magnitudes of the 
effects that are being calculated. But as we will show, 
even the numbers obtained match quite closely with the 
exact numerical results. 

Our idea is to capture the dynamics through a single 
equation ( as far as possible or at most two as in the high 
I and V case) by including all the relevant time scales and 
solve the relevant equation on each branch. For this we 
note that the equation for a and v play a crucial role as 
the inertial contribution appears only through Eqs. © 
and J7J and the time spent by the system is controlled 
by the equation for v, Eq. l(T2|) . Using Eqs. I© and (7J|, 
we get 



F(t)Ra 
a = v/R. 



(15) 



The general equation for a can be written down by using 
Eq. (H2J, in Eq. (JTSJ, we get 



a 



FRa [F(l + a)+ Fa] 

~ W ' 

FRa [F + Fa] 

I Rf 



(16) 
(17) 



In obtaining Eq. (|17fl . we have used 1 + a ~ 1 which 
is valid except for high / and high V. Further, in 
most cases, we can drop Raa as the magnitude of 
this term is small and use F ~ k(V — v). To be 
consistent we use F(t) ~ Fi n + k(V — v)t. We note 
however that even for high / and high V where a is 
not small, dropping l+a and Raa causes only 10% error. 

Case I, small / 

On the low velocity branch AB, as v/R is small in Eq. 
(JSJ), we can drop v term in Eq. 1151 Thus, 



la 



-FRa. 



(18) 



Note that for the low inertia case, sina ~ a approxima- 
tion is clearly justified [see Eq. Q]. Using this equation, 
we first get an idea of the relevant time scales as / is 
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increased. 

Case a 

Consider the low velocity branch AB where the small 
amplitude high frequency oscillations are seen on the 
nearly linearly increasing part of F [given by F(t) = 
F m in + k(V — v)t, see for instance Fig. EJb)]. A rough 
estimate of this time spent on this branch is obtained 
by [fmax - fmin)/kV ~ t. Using f max ~ 284 and 
f mm ~ 200, [from Fig. 13(b)], we get t = 0.084 (com- 
pared to the correct value of 0.063 which we shall ob- 
tain soon) which is much larger than the period of the 
high frequency oscillation. Thus, wc could take the lo- 
cal value F for the purpose of calculating the period 
of the high frequency oscillation. Consider the orbit 
at the lowest value of F for which we can use F m i„ ~ 
/ m i„(i),l) ~ 200. Then using Eq. JTSJ, the frequency 
v = y/(FR/I)/2Tr = 225 for I = 10" 5 which gives the 
period of oscillation T = 4.44 x 10~ 3 . This agrees very 
well with the exact numerical value T = 4.1 x 10~ 3 . This 
frequency decreases when the force reaches the maxi- 
mum value F max ~ f ma x(v, 1) ~ 284 to v = 261 giv- 
ing T — 3.69 x 10 -3 which is again surprisingly close to 
the numerical value 3.72 x 10~ 3 . In the numerical so- 
lutions, we find that the period gradually decreases [see 
Fig. Efc)] . This feature is also easily recovered by using 
F = F m i n + k(V — v)t. This leads to an additional term 
in the equation of motion for a in Eq. (|18fl . 

Ia = -F(t)Ra/I = - [F rnin + k(V - v)t]Ra/I, (19) 

where t is the time required for F to reach F max starting 
from F m i n . Here again the v term can be dropped. If 
Fmin was absent, the equation has the Airy's form. (Note 
that for this case also we could assume F min ~ f m in 
and F max ~ fmax-) Though this equation does not 
have an exact solution, we note that we could take a 
to have a sinusoidal form with 2m/ — yf[FR/I) where 
F is treated as a slowly increasing parameter. (This as- 
sumption works quite well.) The above equation captures 
the essential features of the numerical solution. The nu- 
merical solution of Eq. I|19(l ( as also this representation 
) gives the decreasing trend of the small amplitude high 
frequency oscillations. (Note that the Airy equation itself 
gives a decreasing amplitude 0.) 

We note that Eq. 1(18(1 is valid on the AB branch where 
v is small even for high inertia and small V case. Thus, 
we may be able to recover the gross time scales using 
this equation. Our numerical results show that as we 
increase the inertia, a exhibits a sinusoidal form on the 
AB branch [see Fig. 0(b)], although one full cycle is not 
seen. We note that though the value of a is much larger 
than that for small /, we can still use the above equations 
[Eq. ((18(1 and (|19[l ]. On this branch F increases from a 
value F min ~ fmin to a maximum F max ~ fmax- For 
large I = 10~ 2 ( and V — 1), we get a rough estimate 
of the period by using the mean value of F ~ 240 in Eq. 
PHsl This gives a period T = 0.128 which already agrees 
satisfactorily with the numerically exact value T =0.11 



considering the approximation used (i.e., using the mean 
F). A better estimate can be obtained by using Eq. H19fl . 

For the high / and V case, Fig. 5 for V — 4 shows 
that the wave forms are nearly sinusoidal except for a 
jitter at the top and bottom. For this case, /(w,4) is 
nearly flat over the entire range of values of v, with a 
value ~ 300. Here, even on the AB branch, we can not 
ignore the v term in Eq. 1)15(1. However, one sees that 
as v m i n = 0.335 and v max = 1.25 which suggest that 
to the leading order, we could ignore the v term. This 
gives the period T = 0.115. From Fig. |5fb), considering 
only the monotonically decreasing part (ab) , the value of 
T/2 = 0.53 read off from the figure compares reasonably 
well with this value. 

For the CD branch, as v is not small, the v/R term 
appears to be important in Eq. I|15|) . Some idea of when 
this term is important can be had by looking at the time 
scales arising from inertia, namely, FR/I and the coeffi- 
cient of the damping term, F/Rf in Eq. 1)170. Consider 
V = 1 for I = 10~ 5 and 10~ 2 . The period obtained 
by assuming the mean value of F = 240 in FR/I gives 
4 x 10~ 3 for / = 10~ 5 compared to 0.128 for / = 10~ 2 . 
These numbers can be compared with the time scale 
Rf'/F which is 0.01 (where we have used /' ~ 25 from 
numerical simulations for V = 1). This shows that for 
high inertia the damping coefficient F j Rf in Eq. H17|) is 
important. We will discuss this issue in more detail later. 

Case b 

Now we focus on the origin of jumps between the 
branches. We note that the jumps from CD to AB (or 
vice versa) occur only when the peel velocity v reaches 
a value where f'(v,V) = 0. This also means that the 
time scale on each branch, whether it spends only a short 
time or not, is controlled by the equation for v. However, 
clearly the influence of inertia needs to be included. Here 
we present an approximate equation for v which is valid 
in the various limits of the parameters: 

v = [F(t)(l + a(t))+F(t)a]/f, (20) 
^ [HV ~v) + (Fin + k(V - v)t)a)] ^ 

where the time t is time spent on the branch consid- 
ered ( low or high V). In Eq. 1)21(1. we have again used 
F ~ k(V - v) and F ~ F in + k(V - v)t with the same 
approximation used in Eq. I(17|) . 

We now attempt to obtain correct estimates of the time 
spent by the orbit on each branch starting with the least 
complicated situation of the low inertia and small V. For 
this case, on the low velocity branch, one can use the 
sinusoidal solution for a, namely a = on n sin(2iwt + </>), 
where <f> is a phase factor which also includes the contribu- 
tion arising from the jump as well and 2nv = y/(FR/I) 
with F ~ Fi n + kVt . Both on n and <fi needs to be sup- 
plied. Alternately, one can use Eq. 1(18(1 with Eq. 121() 
for which we provide oti n and cti n at the point from the 
exact numerical solutions. We stress that this procedure 
is not equivalent to solving all the equations, as the only 
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equation we use is Eq. 112111 with the form of a already 
determined from the equation for a. [We note here that 
though we have used the sinusoidal form of a along with 
the initial conditions on ai n ,<j), it is simpler to supply 
the initial conditions cti n , dj„ and use Eq. I|18l) .] We 
note here that /' is a crucial factor that determines the 
time at which the orbit jumps from one branch to the 
other. Equation (|21fl needs to be integrated from Vi n to 
Vf that are determined by the pulling velocity V , i.e., the 
form of f(v,V). 

For the low v branch /' term makes a significant con- 
tribution for the time spent by the trajectory on AB. 
Indeed, one can obtain the order of magnitude of the 
time spent by the orbit on AB by using a crude approx- 
imation for f'(v, 1) = —230(0.5 — v). This can be easily 
integrated from v = Vi n ~ 0.0188, to v = Vf ~ 0.4 which 
already gives At = 0.075. This number is comparable to 
the numerically exact value 0.063. A correct estimate can 
be obtained by using /' from Eq. I|14f> with the sinusoidal 
form a or Eq. I|18|) . (We have used Fj„ = 211.5 from the 
numerical simulations for V = 1 and 2-kv = y/[RF(t)/I] 
with F — Fi n + k(V — v)t.) This gives nearly the exact 
numerical value of At — 0.063. In fact, this solution also 
captures the oscillatory growth nature of v quite accu- 
rately. The approximate form of v(t) (continuous line) 
along with the numerically exact solution ( dotted line) 
are shown in Fig. |5] Using At in F = F in + k(V — v)At 
gives AF = 63 and F = 274.5 which is in good agree- 
ment with the exact numerical value of F max = 275. It 
is interesting to note that this value is much less than 
fmax — 283 [see Fig. Etb)] or equivalently AF is less than 
fmax ~ /mim what is also observed in our exact numeri- 
cal simulation. The underlying mechanism of jumping of 
the orbit before F reaches fmax also becomes clear from 
the analysis (Fig. [5J. We note that the magnitude of the 
oscillatory component in v grows till it reaches v max per- 
mitted by f(v, 1). Then, the orbit has to jump to CD. 
Thus, the approximate solution gives an insight into the 
cause of the orbit jumping even before F reaches fmax 
(for small I). 

For the CD branch also, the dominant term is /'. In- 
deed, any reasonable function which has the same geo- 
metrical form of / shown in Fig. [5] will give good results 
for At. Using the correct form of /', we get At = 0.005 
which is close to the exact result. This again gives cor- 
rect magnitude of AF = 72.5. In addition the nature of 
the v(t) obtained by this approximation is close to the 
exact numerical solution shown in the inset of Fig. 8. 
Case II, intermediate and high I and low V 

The most difficult feature of our numerical solutions 
to understand is the dynamical mechanism leading to a 
series of drops in the pull force seen on the descending 
branch of F(t) for intermediate and high values of inertia 
and for a range of V values. Consider the high inertia and 
low V case (say I — 10~ 2 and V = 1) shown in Fig. 0] 
As stated earlier, there are two different issues that need 
to be understood here. First, the series of small force 
drops AF and second the monotonic increasing nature 
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FIG. 8: Comparison of approximate solution (continuous) 
with the numerically exact solution (dotted) of v(t) for I — 
10 -5 and V — 1 for the AB branch. The inset shows a similar 
comparison of v(t) on the CD branch, [v, V are in m/s, / in 
kg m 2 and t in s.) 



of F on the AB branch. 

In this case, as already discussed, the coefficient of a, 
namely, F/Rf term in Eq. (|17J) determines the time 
scale on CD, while on AB, the term FRa/I dominates. 
Thus, the general equation valid for this case is 



FRa 



Fa 



(22) 



where we use F = Fi n + k(V — v)t. [Note that we have 
dropped k(V — v) term from Eq. (|17f) as this term does 
not have any dependence on a or a.] 

We start with the cascading effect. Consider the orbit 
when it is at the highest value of Fj„ = 295.6 on the CD 
branch for which we can drop FRa/I term. As /' is a 
function of v, and F also depends on time, it appears that 
we need to use coupled equations a = —Fa/Rf with Eq. 
I)2ip. However, the numerical solution of these equations 
show that one can make further approximation by taking 
/' to be constant taken at v = 15.54 and F = Fj„, as 
the time spent on this branch is very small. The error 
in using this approximation is within 10%. Indeed, using 
ai n = —0.0304, ctin = —160 and numerically integrat- 
ing Eq. J2U, along with Eq. from v ln = 15.54 to 
Vf = 8.7 gives At = 1.78 x 10~ 3 . This compares reason- 
ably with the numerical value of 1.89 x 10~ 3 . Using this 
we get AF = 19.4 which compares well with the numer- 
ically exact value 19.8. At this point the orbit jumps to 
the low velocity AB branch (to the point e). Thus, as 
At is small, for all practical purposes, we can ignore the 
dependence of / on v and F on t and use a to be an ex- 
ponentially decreasing function for analytical estimates. 
These analytical estimates already give reasonably accu- 
rate numbers. 

On the AB branch, the dominant time scale is deter- 
mined by FRa/I, and we can use the approximate sinu- 
soidal form in Eq. lj2T|) , or Eq. ftjgjl along with Eq. J2U 



11 



for the time evolution from the point e. Integrating from 
v in = 0.0188 to Vf = 0.4 with the appropriate initial val- 
ues ai n = 0.239, a i n = 11.9 (or a, 4>) and F in — 276.53, 
gives At = 0.016 which again compares very well with ex- 
act numerical value At — 0.0164. This gives AF = 11.61. 
The procedure for calculating the time spent by the orbit 
on CD and AB is the same and we find that successive 
values of AF increases which is again consistent with 
what is seen in Figs. Ufa) andQJc). 

Continuing this procedure, we find that a minimum 
value of F — 186.95 for the cycle is reached. Now 
consider the time evolution of F on AB that should 
lead to a monotonically increasing nature as seen in the 
numerically exact solution. As this point corresponds 
to the point at which the dynamics switches from the 
jumping mode to the monotonically increasing nature 
of F (i.e., the stretch ab), we discuss this in some de- 
tail. For the point a, we have used the initial condition 
a in — 0.0599, ai n = 9.7 and integrating Eq. l|2*T|) and Eq. 
(jTBjl (or the sinusoidal form of a) from v = v min = 0.0188 
to v = v max = 0.4 gives At — 0.117. This is nearly the 
value 0.114 obtained from the exact numerical integra- 
tion. This gives AF = 117 and F ma x = 303.95 which 
compares very well with the exact numerical value. In 
addition, the growth form of v obtained from this ap- 
proximation (continuous line) agrees very well with that 
of the exact numerical solution (dotted line) as shown in 
Fig. [5J The discrepancy seen in the figure can be reduced 
for instance if we include the terms neglected in Eq. H16|) 
such as F and using 1 + a in Eq. (|2*U)l . 

Now we come to the crucial question. How does the 
system know that it has to go from a to b, while just 
during the previous visit to the point k on AB branch 
lead only to a small increase in AF [Fig. 0Ja) and^fc)] 
before jumping to CD? To understand this, we recall 
that on AB, a sinusoidal solution is allowed. First, one 
can notice a few differences in the initial conditions 
between the point a and k. For the point k, the initial 
conditions taken from the exact numerical solution are 
a in = 0.298 and a in = 18.3. ( F in = 193.37), while for 
the point a, a,„ = 0.0589, a,-„ = 9.7. However, for a 
to begin a sinusoidal form, the initial value of a = 18.3 
is much higher than the natural slope. The local slope 
for any sinusoidal form is maximum when the variable 
is close to zero. In Fig. 0Jb), the sinusoidal form starts 
when a is close to zero (~ 0.0589 at a) where the local 
slope should be close to the maximum value. Near 
a ~ 0, the local slope is the product of the maximum 
amplitude of a, say, cto (in the sinusoidal stretch ab) and 
2m/. (We have assumed a = a$ sinl-nvt by dropping the 
phase factor.) Thus, one should have a ~ l-nva^ when 
a ~ 0. Using the value ao = 0.23 from exact numerical 
solution and v ~ 6.88 at Fj„ = 186.95, we find that 
a ~ 10 near a ~ 0. Indeed, this is satisfied only at a 
where 6n n — 9.7. [Note that a is not symmetric around 
zero due to the presence of v in Eq. © which has been 
ignored for the purpose of present discussion.] However, 
ain — 18.3 at k is significantly higher than the slope 



permitted for a to start a sinusoidal sojourn. This forces 
the orbit to make one more small loop (AB to CD and 
back) so that the initial value of &i n is commensurate for 
a to start a sinusoidal form. Indeed, the initial values of 
a at all the earlier visits to AB branch keep decreasing 
until it reaches a value that is consistent to begin the 
sinusoidal growth. Once this is satisfied, the monotonic 
increasing behavior from a to 6 is seen. As we will show 
this is the mechanism operating for high I and V case. 

Case III, High I and V 

For this case, even on the AB branch, v/R cannot be 
ignored in Eq. (|T5)l and thus one needs to use coupled 
Eqs. (|21|l and (|15l) . Calculations follow much the same 
lines and give correct values for At and AF on both the 
branches during the rapid jumps. 

Again, we need to answer when exactly does the system 
know to switch from a rapid jumping mode to monoton- 
ically increasing on AB or decreasing mode on CD? 

Consider the last of the rapid jumps from CD to AB 
(just prior to the point a) in Fig. |5fc). The corresponding 
point in the a plot [Fig. E£b)] is shown on an expanded 
scale in the inset. From this figure, it is clear that a has 
a positive slope at k, though of small magnitude while 
at a, it has a value -9.7. The latter is close to the nat- 
ural (negative) slope of a when it begins the descending 
branch of the sinusoidal form. On the other hand, the 
slope of a is positive at k and hence will not allow the 
growth to change over from a jumping mode to the sinu- 
soidal growth form for a. One can note that the slopes 
at points of all the earlier visits to AB [see Fig. E^b) 
inset] keep decreasing till the slope becomes negative re- 
quired for the monotonically decreasing trend of a. This 
is exactly the same mechanism for I = 10~ 2 and V = 1 
also, for the low v branch, except that in this case, even 
the sign of the slope is incompatible for all the points 
prior to a in Fig. [5]c. The mechanism operating on CD 
(i.e., at the switching from jumping mode to monotoni- 
cally decreasing nature of F) is essentially the same but 
arguments are a little more involved and hence they are 
not presented. 

Now, we consider the causes leading to the maximum 
and minimum values taken by F being much more than 
permitted by f(v,V). As this is dominant for V = 4, 
we illustrate this using Fig. |3fa) and 0c). We first note 
that Eq. HI U| l constrains the dynamically changing val- 
ues of F(t) and ait) to the stationary values of f(v, V). 
Clearly, this implies that F = f(v, V)/(l + sina). A 
rough estimate of F max can be obtained by F max ~ 
fmax/ (1+Omm) with F mm determined by a max - This re- 
lation can be easily verified by using the numerical values 
of a. For instance, for V — 4 and I = 10~ 2 , a m in = —0.3 
and fmax = 307. This gives F m ax — 438 while the nu- 
merical value from the phase plot for this case gives 433 
which is very close. Similarly, using a ma x = 0.62 and 
fmin = 293, we get F m i n = 181 which compares well 
with the numerical value of 180. We have verified this 
relation is respected for various values of V and I. For 
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FIG. 9: Comparison of the approximate solution (continuous 
line) for v(t) with numerically exact solution (dotted line) for 
V = 1,1 = 10 -2 . (v, V are in m/s, I in kg m 2 and t in s.) 



small I, a is small, we should not find much difference 
between F max (F min ) and that of /. 



VI. SUMMARY AND CONCLUSIONS 

We first summarize the results before making some 
relevant remarks. We have carried out a study of the 
dynamics of an adhesive roller tape using a differential- 
algebraic scheme used for singular set of differential equa- 
tions. The algorithm produces stick-slip jumps across the 
two dissipative branches as a consequence of the inher- 
ent dynamics. Our extensive simulations show that the 
dynamics is much richer than anticipated earlier. In par- 
ticular the influence of inertia is shown to be dramatic. 
For instance, even at low inertia, for small values of V, 
the influence of inertia manifests with jumps of the orbit 
occurring even before F reaches f m ax (or f m in) which 
is quite unexpected. More dominant is its influence for 
high I both for low V and high V, though it is striking 
for the latter case. Following the reasoning used in the 
PLC effect, we introduce a dynamized curve f(v, V) as 
resulting from competing time scales of internal relax- 
ation and imposed pull speed. The modified peel force 
function leads to the decreasing trend in the magnitude 
of AF with increasing pull velocity, a feature observed in 
experiments. We have also recaptured the essential fea- 
tures of the dynamics by a set of approximations valid in 
different regimes of the parameter space. These approx- 
imate solutions illustrate the influence of various time 
scales such as that due to inertia, the elasticity of the 
tape and that determined by the stationary peel force 
f(v,V). We also find the unusual canard type of solu- 
tions. 

Here, it is worthwhile to comment on the dynamical 
features of the model. The numerical results themselves 
are too complex to understand. A striking example of 
this is the series of force drops seen on the descending 
branch of the pull force [Fig. Etc)]. This result is hard 
to understand as it would amount to a partial relaxation 



of the pull force. However, a partial relaxation is only 
possible in the presence of another competing time scale 
(other than the imposed time scale). Another example is 
the jumping of the orbit for low I case, from AB to CB 
and vice versa even before the pull force reaches the ex- 
tremum values of f[v, V). For this reason, we have under- 
taken to make this complex dynamics transparent using 
a set of approximations. The basic idea here is to solve a 
single equation (or at most two equations as in the high 
/ and V case) which incorporates all the relevant time 
scales. This method not only captures all the results to 
within 10% error but it also clearly brings out the regimes 
of parameter space where these time scales become im- 
portant. This analysis also shows that the time scale due 
to inertia of the roller tape shows up even for low I which 
comes as a surprise as one expects that for low inertia, 
the orbit should stick to the stationary peel function. 
(Recall that for low inertia, equations have been approx- 
imated by Lienard type of equations by Maugis and Bar- 
quins 0.) Our approximate equations demonstrate that 
a crucial role in inducing the jumps even at low inertia is 
played by the high frequency oscillations resulting from 
the inertia of the roller tape. As for high inertia (both 
for low and high pull velocities), the time scale due to 
inertia is responsible for the partial relaxation of F as 
shown. 

A few comments may be in order on the bursting type 
of oscillations in the peel velocity. Bursting type of os- 
cillatory behavior are commonly seen in neuro-biological 
systems [25|. Conventionally, bursting typ e oscillations 
arise in the presence of homoclinic orbit [25| . Such burst- 
ing type of oscillations have also been modeled using one 
dimensional map [26j. However, it is clear that the mech- 
anism for bursting type of oscillations in our case is dif- 
ferent. In our case, this arises due to the fact that the 
orbit is forced to jump between the stable manifolds as 
a result of competing time scale of inertia and the time 
scale for the evolution of v. (We note that the latter itself 
includes more than one time scale [see Eq. (I21H ]. namely 
the contribution from the slopes of the stable parts of the 
stationary curve f(v, V) and that due to elasticity of the 
tape. ) The bunching of the spikes in v is the result of 
f{v,V) becoming flat for large V and /. One other com- 
ment relates to canard type solutions. Figure \J\ shows 
one such solution. As mentioned, these type solutions 
arise from sticking to the unstable manifold. In fact, a 
similar type of solution is seen in Fig. |5fc). As noted ear- 
lier, all the jumps from CD to AB or vice versa always 
occur when the peel velocity reaches the limiting value 
where f'(v,V) = 0. However, it can be seen from this 
figure, the orbit starting from c monotonically decreases 
well into the unstable part of /. Thus, this solution also 
has the features of canards. It must be stated that our 
approximate solutions cannot capture the behavior of ca- 
nards. 

Finally, the results presented in this paper are on the 
nature of dynamics of the model equations which so far 
had defied solution. However, comparison with exper- 
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iments has been minimal largely due to the paucity of 
quantitative experimental findings as stated earlier. Our 
analysis shows that the model predicts periodic, saw- 
tooth as well as chaotic solutions as reported in 
The high magnitude of the Lyapunov exponents for the 
chaotic solutions in the lowpull velocities is consistent 
with that reported earlier jj]. We note that the other 
quantitative experimental feature reported by Refs. 0, Q 
is the decreasing trend of the average force drop magni- 
tudes as a function of the pull velocity is also captured by 
our model (Fig. EJ), a result that holds for both low and 
high inertia. This result is a direct consequence of the 
dynamization of the peel force function, i.e., dependence 
of the peel force on the pull velocity. We note here that 
the complex dynamics at high velocities (see Fig. EJl is 
a direct result of the unstable part of dynamized curve, 
f(v, V), shrinking to zero. To the best of our knowledge, 
this is first time the result in Fig. [S]has been explained. 
As the hypothesis of dynamization captures the decreas- 
ing trend of the force drops, it also suggests that the 



underlying mechanism of competing time scales respon- 
sible for the peel force depending on the pull velocity is 
likely to be correct as in the PLC effect. Clearly, a rigor- 
ous derivation of the peel force function from microscopic 
considerations that includes the effect of the viscoelastic 
glue at the contact point is needed to understand the 
dynamics appropriately. 
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